Fungal OTU’s across life stages of boreal toad: published and supplemental analyses
Alexandra Alexiev
12/17/2019
#packages and set working directory
library(vegan)
library(ggplot2)
library(mctoolsr)
library(tidyr)
library(tibble)
library(gridExtra)
library(dplyr)
library(ggpubr)
library(cowplot)
library(kableExtra)
library(RColorBrewer)
set.seed(33333)
work_dir <- "/Users/alal3493/Documents/Projects/DevelopmentalBorealToad/FungalDevData/02_PublishedAnalysis/"
seq_dir <- "/Users/alal3493/Documents/Projects/DevelopmentalBorealToad/FungalDevData/00_SeqProcessing/"Basic measures of highly abundant or present taxa
#load in starting files
mapfp_all <- paste0(work_dir,"Allsamps_clean_maptab.txt")
tabfp_all <- paste0(work_dir,"Allsamps_clean_otutab.txt")
input_all <- load_taxa_table(tab_fp = tabfp_all, map_fp = mapfp_all) # this file has both toad and environment sample types in it (169 samples)
mapfp_t <- paste0(work_dir,"Toad_clean_maptab.txt")
tabfp_t <- paste0(work_dir,"Toad_clean_otutab.txt")
input_toad <- load_taxa_table(tab_fp = tabfp_t, map_fp = mapfp_t) # this file only has toad samples in it (124 samples)
# order by life stage
life_list <- c("Sediment", "Water",
"Eggs.11.15", "Tadpole.20.22",
"Tadpole.23.25", "Tadpole.25.27",
"Tadpole.27.29", "Tadpole.29.31",
"Tadpole.31.35", "Tadpole.36.39",
"Metamorph.40.46", "Subadult", "Adult")
input_all$map_loaded$Gosner_Stage <-
factor(input_all$map_loaded$Gosner_Stage,
levels = life_list)
input_toad$map_loaded$Gosner_Stage <-
factor(input_toad$map_loaded$Gosner_Stage,
levels = life_list)
# and order the other column I'll be using as well, where tadpoles are lumped
samptype_list <- c("Sediment", "Water", "Eggs",
"Tadpole", "Metamorph", "Subadult", "Adult")
input_all$map_loaded$Life_Stage_Simplified <-
factor(input_all$map_loaded$Life_Stage_Simplified,
levels = samptype_list)
input_toad$map_loaded$Life_Stage_Simplified <-
factor(input_toad$map_loaded$Life_Stage_Simplified,
levels = samptype_list)# number of total fungal taxa found in all samples
nrow(input_all$data_loaded) # 2156
## [1] 2156
nrow(input_toad$data_loaded) # 1760
## [1] 1760
# top ten taxa on toads and in environment
toptaxa_all <- return_top_taxa(input = input_all, number_taxa = 10)
toptaxa_toad <- return_top_taxa(input = input_toad, number_taxa = 10)
kable(toptaxa_all, caption = "Top ten taxa in toad and environment samples") %>%
kable_styling(bootstrap_options = c("striped","bordered"),
full_width = F, position = "left")| taxonomy1 | taxonomy2 | taxonomy3 | taxonomy4 | taxonomy5 | taxonomy6 | taxonomy7 | |
|---|---|---|---|---|---|---|---|
| OTU_1 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Capnodiales | f__Mycosphaerellaceae | g__Mycosphaerella | s__Mycosphaerella_tassiana |
| OTU_2 | k__Fungi | NA | NA | NA | NA | NA | NA |
| OTU_4 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Capnodiales | f__Cladosporiaceae | g__Cladosporium | s__Cladosporium_delicatulum |
| OTU_3 | k__Fungi | p__Ascomycota | c__Leotiomycetes | o__Helotiales | f__Hyaloscyphaceae | g__Pezizella | s__Pezizella_epithallina |
| OTU_7180 | k__Fungi | p__Basidiomycota | c__Agaricomycetes | o__Agaricales | f__Schizophyllaceae | g__Schizophyllum | s__Schizophyllum_commune |
| OTU_8 | k__Fungi | p__Ascomycota | c__Leotiomycetes | o__Thelebolales | f__Pseudeurotiaceae | g__Pseudeurotium | s__Pseudeurotium_hygrophilum |
| OTU_34 | k__Fungi | p__Rozellomycota | c__unidentified | o__unidentified | f__unidentified | g__unidentified | s__unidentified |
| OTU_11 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Pleosporales | f__Phaeosphaeriaceae | g__unidentified | s__unidentified |
| OTU_12 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Pleosporales | f__Phaeosphaeriaceae | g__Sclerostagonospora | s__Sclerostagonospora_phragmiticola |
| OTU_37 | k__Fungi | p__Basidiomycota | c__Tremellomycetes | o__Filobasidiales | f__Filobasidiaceae | g__Naganishia | s__Naganishia_albida |
kable(toptaxa_toad, caption = "Top ten taxa in toad samples") %>%
kable_styling(bootstrap_options = c("striped","bordered"),
full_width = F, position = "left")| taxonomy1 | taxonomy2 | taxonomy3 | taxonomy4 | taxonomy5 | taxonomy6 | taxonomy7 | |
|---|---|---|---|---|---|---|---|
| OTU_1 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Capnodiales | f__Mycosphaerellaceae | g__Mycosphaerella | s__Mycosphaerella_tassiana |
| OTU_4 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Capnodiales | f__Cladosporiaceae | g__Cladosporium | s__Cladosporium_delicatulum |
| OTU_2 | k__Fungi | NA | NA | NA | NA | NA | NA |
| OTU_7180 | k__Fungi | p__Basidiomycota | c__Agaricomycetes | o__Agaricales | f__Schizophyllaceae | g__Schizophyllum | s__Schizophyllum_commune |
| OTU_11 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Pleosporales | f__Phaeosphaeriaceae | g__unidentified | s__unidentified |
| OTU_3 | k__Fungi | p__Ascomycota | c__Leotiomycetes | o__Helotiales | f__Hyaloscyphaceae | g__Pezizella | s__Pezizella_epithallina |
| OTU_15 | k__Fungi | p__Basidiomycota | c__Tremellomycetes | o__Filobasidiales | f__unidentified | g__unidentified | s__unidentified |
| OTU_34 | k__Fungi | p__Rozellomycota | c__unidentified | o__unidentified | f__unidentified | g__unidentified | s__unidentified |
| OTU_37 | k__Fungi | p__Basidiomycota | c__Tremellomycetes | o__Filobasidiales | f__Filobasidiaceae | g__Naganishia | s__Naganishia_albida |
| OTU_12 | k__Fungi | p__Ascomycota | c__Dothideomycetes | o__Pleosporales | f__Phaeosphaeriaceae | g__Sclerostagonospora | s__Sclerostagonospora_phragmiticola |
Q1: What are the factors associated with variation in fungal communities?
Alpha diversity differences between life stages and environment samples
Below is a supplementary figure comparing three alpha diversity measures across all the sample types. Since site is a confounding factor here, given exploratory analyses, I analyzed with blocking by site (“strata”).
# calculate alpha diversity measures
input_all$map_loaded$rich = calc_diversity(
tax_table = input_all$data_loaded, metric = 'richness')
input_all$map_loaded$simpson = calc_diversity(
tax_table = input_all$data_loaded, metric = 'simpson')
input_all$map_loaded$shannon = calc_diversity(
tax_table = input_all$data_loaded, metric = 'shannon')
# Kruskall-Wallis test on life stages, strata-site
# species richness
krusk_rich <- kruskal.test(input_all$map_loaded$rich ~
input_all$map_loaded$Gosner_Stage)
bonpval_rich <- p.adjust(p = krusk_rich$p.value, method = "bonferroni",
n = length(unique(input_all$map_loaded$Gosner_Stage)))
krusk_rich; bonpval_rich
##
## Kruskal-Wallis rank sum test
##
## data: input_all$map_loaded$rich by input_all$map_loaded$Gosner_Stage
## Kruskal-Wallis chi-squared = 102.79, df = 12, p-value < 2.2e-16
## [1] 2.048359e-15
# simpson (evenness)
krusk_simp <- kruskal.test(input_all$map_loaded$simpson ~
input_all$map_loaded$Gosner_Stage)
bonpval_simp <- p.adjust(p = krusk_simp$p.value, method = "bonferroni",
n = length(unique(input_all$map_loaded$Gosner_Stage)))
krusk_simp; bonpval_simp
##
## Kruskal-Wallis rank sum test
##
## data: input_all$map_loaded$simpson by input_all$map_loaded$Gosner_Stage
## Kruskal-Wallis chi-squared = 29.706, df = 12, p-value = 0.003092
## [1] 0.04019437
# shannon (both richness and evenness)
krusk_shan <- kruskal.test(input_all$map_loaded$shannon ~
input_all$map_loaded$Gosner_Stage)
bonpval_shan <- p.adjust(p = krusk_shan$p.value, method = "bonferroni",
n = length(unique(input_all$map_loaded$Gosner_Stage)))
krusk_shan; bonpval_shan
##
## Kruskal-Wallis rank sum test
##
## data: input_all$map_loaded$shannon by input_all$map_loaded$Gosner_Stage
## Kruskal-Wallis chi-squared = 52.188, df = 12, p-value = 5.741e-07
## [1] 7.463869e-06
# GRAPHING
# Species richness graphs: mean number of OTUs found per sample type
liferich <- ggplot(input_all$map_loaded, aes(Gosner_Stage, rich)) +
geom_boxplot() + xlab("") +
ylab("Fungal OTU Richness") + ggtitle("A") +
geom_text(x = 7, y = 194, size = 3,
label = paste0("KW chi-sq stat=",signif(krusk_rich$statistic,
digits = 5),
"\n","p=",signif(bonpval_rich, digits = 3))) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, colour = "black"),
axis.text.y = element_text(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"))
# Simpson's: measures species/OTU evenness, or how close in number each species in the environment
lifesimp <- ggplot(input_all$map_loaded, aes(Gosner_Stage, simpson)) +
geom_boxplot() + xlab("Sample Type") +
ylab("Simpson's Diversity Metric") + ylim(0,1.2) + ggtitle("B") +
geom_text(x = 7, y = 1.15, size = 3,
label = paste0("KW chi-sq stat=",signif(krusk_simp$statistic,
digits = 5),
"\n","p=",signif(bonpval_simp, digits = 3))) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, colour = "black"),
axis.text.y = element_text(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"))
# Shannon: combination of richness and evenness
lifeshan <- ggplot(input_all$map_loaded, aes(Gosner_Stage, shannon)) +
geom_boxplot() + xlab("") +
ylab("Shannon Diversity Index") + ggtitle("C") +
geom_text(x = 7, y = 4.1, size = 3,
label = paste0("KW chi-sq stat=",signif(krusk_shan$statistic,
digits = 5),
"\n","p=",signif(bonpval_shan, digits = 3))) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, colour = "black"),
axis.text.y = element_text(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"))
alpha_panel <- grid.arrange(liferich, lifesimp, lifeshan, ncol = 3)
alpha_panel
## TableGrob (1 x 3) "arrange": 3 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (1-1,3-3) arrange gtable[layout]
# ggsave(paste0(work_dir,"Figs_Tables/AllAlphaDiv.png"),
# alpha_panel, dpi = 300,
# width = 10, height = 7)Pairwise OTU richness Stats
- Here, I’m using a two-sample Mann-Whitney test (or unpaired wilcoxon test). This is a non-parametric test of whether two populations are different across two different life stages, while still assuming those populations are independent. Since my samples came from different animals at different life stages, I need to make this assumption, so I cannot use a paired wilcoxon test.
- I calculated Bonferroni (signif if p<=0.05) and FDR (p<=0.05) p-value corrections, used both together to figure out how big of a difference there was between each life stage.
# Are pairwise alpha diversities same or different across life stages?
nlevels(input_all$map_loaded$Gosner_Stage) # 13
## [1] 13
cats <- levels(input_all$map_loaded$Gosner_Stage) # define category levels
cats_pairs <- combn(cats, 2) # make pairs of these
pvals <- c() # output df for pvalues
Wstat <- c() # output df for Wstat
for ( i in 1:ncol(cats_pairs) ) { # for each pair of sample types
pair <- cats_pairs[, i] # define the pair
Gosner_rich <- input_all$map_loaded %>% # take the mapping file
dplyr::select(Gosner_Stage, rich) %>% # select these two columns from it
dplyr::filter(Gosner_Stage %in% pair) # filter the sample types in the pair
wil <- wilcox.test(rich ~ Gosner_Stage, data = Gosner_rich, paired = F, exact = F) # run Wilcoxon test
pvals <- c(pvals, wil$p.value) # save the p val
Wstat <- c(Wstat, wil$statistic) # save the W statistic
}
results <- data.frame(t(cats_pairs), Wstat, pvals) # make a dataframe of the useful outputs
results$pvalBon = pvals * length(pvals) # calc bonferroni pval correction
results$pvalFDR = round(pvals * (length(pvals) / rank(pvals, ties.method = "average")),
3) # calc FDR pval correction
results <- results %>%
arrange(pvalBon, pvalFDR) %>% # arrange in ascending order by the corrected pvals
filter(pvalBon <= 0.05) %>%
filter(pvalFDR <= 0.05) %>% # filter by the two corrected pvals being <= 0.05
dplyr::select(-pvals) # remove pval column
# write.table(results, paste0(work_dir, "Figs_Tables/pairedMW_alphadiv.txt"),
# sep="\t", quote = F, row.names = F)Make three paneled figure with alpha and beta diversity measures
# make color palette
nlevels(input_all$map_loaded$Life_Stage_Simplified)
## [1] 7
brewer.pal(n = 7, name = "Dark2") # display hexcodes for the pallette I'm using but want to change slightly
## [1] "#1B9E77" "#D95F02" "#7570B3" "#E7298A" "#66A61E" "#E6AB02" "#A6761D"
samptype_cols <- c("#A6761D",
"#3a6b94", "#66A61E",
"#E7298A", "#7570B3",
"#D95F02", "#1B9E77")
# make OTU richness figure
OTUrich <- ggplot(input_all$map_loaded, aes(Gosner_Stage, rich)) +
geom_boxplot(alpha = 0.6, outlier.shape = NA,
aes(fill = Life_Stage_Simplified)) + xlab("") +
ylab("Fungal OTU Richness") + ggtitle("A") +
geom_point(aes(fill = Life_Stage_Simplified), size = 3, shape = 21,
position = position_jitter(
width = 0.25
)) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, colour = "black"),
axis.text.y = element_text(colour = "black"),
axis.title = element_text(size = 16),
axis.text = element_text(size = 14),
legend.title = element_text(size = 16),
legend.text = element_text(size = 14),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
plot.title = element_text(size = 30),
legend.position = "none") +
labs(fill = "Sample Type") +
scale_fill_manual(values = samptype_cols) +
scale_x_discrete(labels = c("Sediment", "Water",
"Eggs", "Tadpole (20-22)",
"Tadpole (23-25)", "Tadpole (25-27)",
"Tadpole (27-29)", "Tadpole (29-31)",
"Tadpole (31-35)", "Tadpole (36-39)",
"Metamorph (40-46)", "Subadult", "Adult"))
# make NMDS plot with all samples, toad vs. env
# create distance matrix
md_transformed <- t(sqrt(input_all$data_loaded))
dm <- vegdist(md_transformed, method = "bray")
# run NMDS
md.nmds <- metaMDS(dm, k = 3, trymax = 1000) #solution reached, stress=0.195
## Run 0 stress 0.195206
## Run 1 stress 0.1985572
## Run 2 stress 0.196332
## Run 3 stress 0.195203
## ... New best solution
## ... Procrustes: rmse 0.002155213 max resid 0.02202525
## Run 4 stress 0.1956287
## ... Procrustes: rmse 0.009492756 max resid 0.09631747
## Run 5 stress 0.1969824
## Run 6 stress 0.1954464
## ... Procrustes: rmse 0.01550547 max resid 0.1772666
## Run 7 stress 0.1963652
## Run 8 stress 0.1964445
## Run 9 stress 0.1993956
## Run 10 stress 0.2113578
## Run 11 stress 0.1961069
## Run 12 stress 0.1963005
## Run 13 stress 0.1974691
## Run 14 stress 0.195507
## ... Procrustes: rmse 0.01334855 max resid 0.1575973
## Run 15 stress 0.1964108
## Run 16 stress 0.1952052
## ... Procrustes: rmse 0.0008926742 max resid 0.005093944
## ... Similar to previous best
## Run 17 stress 0.1984376
## Run 18 stress 0.1974831
## Run 19 stress 0.1962913
## Run 20 stress 0.1957653
## *** Solution reached
#prepare NMDS points and values for graphing
md.nmds.points <- md.nmds$points %>% # take the NMDS points
data.frame() %>% # make a data frame
rownames_to_column(var = "mysamples") %>%
mutate(SampleID = mysamples) %>%
column_to_rownames(var = "mysamples")
meta <- input_all$map_loaded %>%
rownames_to_column("SampleID")
md.nmds.metadata <- inner_join(x = md.nmds.points, y = meta,
by = "SampleID") %>%
group_by(Type) %>%
ungroup() %>%
dplyr::select(SampleID, MDS1, MDS2, everything()) # select these columns only
md.nmds.metadata.unq <- md.nmds.metadata %>%
dplyr::select(Type) %>%
unique()
md.nmds.metadata.unq$Type <-
factor(md.nmds.metadata.unq$Type,
levels = c("Toad", "Environment"))
# get stress value
stress.md = paste("stress =", round(md.nmds$stress, digits = 4))
# get NMDS hulls
group.chulls <- plyr::ddply(md.nmds.metadata, "Type",
function(df) df[chull(df$MDS1, df$MDS2), ])
# make color pallete
samptype_cols_TE <- c("#4E72E2", "#E2A14E")
# Plot NMDS
nmds.plot <- ggplot(md.nmds.metadata, aes(x = MDS1, y = MDS2, color = Type)) +
geom_polygon(data = group.chulls,
aes(fill = Type),
alpha = 0.15, linetype = 0) +
geom_point(size = 3) +
annotate("text", x = Inf, y = Inf, label = stress.md, hjust = 1, vjust = 1) +
ggtitle("C") +
labs(colour = "Sample Type") +
scale_colour_manual(values = samptype_cols_TE) +
scale_fill_manual(values = samptype_cols_TE) +
guides(fill = F, color = guide_legend(order = 1,
override.aes = list(shape = 15,
size = 6))) +
theme(axis.title = element_text(size = 16),
axis.text = element_text(size = 14),
legend.title = element_text(size = 16),
legend.text = element_text(size = 14),
plot.title = element_text(size = 30))
# make another NMDS with only the toad life stages shown
# create distance matrix
md_transformed <- t(sqrt(input_toad$data_loaded))
dm <- vegdist(md_transformed, method = "bray")
# run NMDS
md.nmds <- metaMDS(dm, k = 3, trymax = 1000) #solution reached, stress=0.198
## Run 0 stress 0.1990207
## Run 1 stress 0.1997202
## Run 2 stress 0.1994624
## ... Procrustes: rmse 0.02582633 max resid 0.1347371
## Run 3 stress 0.2023672
## Run 4 stress 0.2006704
## Run 5 stress 0.2008573
## Run 6 stress 0.1994031
## ... Procrustes: rmse 0.02100356 max resid 0.120204
## Run 7 stress 0.2001523
## Run 8 stress 0.2051994
## Run 9 stress 0.1987084
## ... New best solution
## ... Procrustes: rmse 0.01775579 max resid 0.1122644
## Run 10 stress 0.1992027
## ... Procrustes: rmse 0.03218216 max resid 0.1800766
## Run 11 stress 0.2020069
## Run 12 stress 0.2004879
## Run 13 stress 0.1990802
## ... Procrustes: rmse 0.01716068 max resid 0.1245285
## Run 14 stress 0.2018043
## Run 15 stress 0.2050105
## Run 16 stress 0.1991984
## ... Procrustes: rmse 0.02372536 max resid 0.1858185
## Run 17 stress 0.1997273
## Run 18 stress 0.2002501
## Run 19 stress 0.205055
## Run 20 stress 0.1990654
## ... Procrustes: rmse 0.02649491 max resid 0.1637236
## Run 21 stress 0.1985955
## ... New best solution
## ... Procrustes: rmse 0.0137297 max resid 0.1201587
## Run 22 stress 0.1987048
## ... Procrustes: rmse 0.009884085 max resid 0.05074342
## Run 23 stress 0.2002589
## Run 24 stress 0.2046748
## Run 25 stress 0.2042179
## Run 26 stress 0.1994072
## Run 27 stress 0.1987614
## ... Procrustes: rmse 0.02432677 max resid 0.1619859
## Run 28 stress 0.2037788
## Run 29 stress 0.1985588
## ... New best solution
## ... Procrustes: rmse 0.007976802 max resid 0.04539034
## Run 30 stress 0.2019626
## Run 31 stress 0.1988353
## ... Procrustes: rmse 0.0262126 max resid 0.1735216
## Run 32 stress 0.2029194
## Run 33 stress 0.1991735
## Run 34 stress 0.2036507
## Run 35 stress 0.2002601
## Run 36 stress 0.2020978
## Run 37 stress 0.1997558
## Run 38 stress 0.1986237
## ... Procrustes: rmse 0.02431091 max resid 0.1718288
## Run 39 stress 0.2050813
## Run 40 stress 0.2005529
## Run 41 stress 0.1990195
## ... Procrustes: rmse 0.02043965 max resid 0.1854972
## Run 42 stress 0.2044782
## Run 43 stress 0.202236
## Run 44 stress 0.199188
## Run 45 stress 0.1993085
## Run 46 stress 0.2013434
## Run 47 stress 0.1990424
## ... Procrustes: rmse 0.02653854 max resid 0.1790137
## Run 48 stress 0.2023573
## Run 49 stress 0.1994305
## Run 50 stress 0.1998933
## Run 51 stress 0.1988918
## ... Procrustes: rmse 0.01179392 max resid 0.05615324
## Run 52 stress 0.2003017
## Run 53 stress 0.1991179
## Run 54 stress 0.2042435
## Run 55 stress 0.1985242
## ... New best solution
## ... Procrustes: rmse 0.007902692 max resid 0.06864319
## Run 56 stress 0.1991631
## Run 57 stress 0.2039104
## Run 58 stress 0.2049592
## Run 59 stress 0.2007415
## Run 60 stress 0.2023174
## Run 61 stress 0.2059404
## Run 62 stress 0.2036642
## Run 63 stress 0.2020151
## Run 64 stress 0.1992684
## Run 65 stress 0.2011889
## Run 66 stress 0.2010427
## Run 67 stress 0.1998556
## Run 68 stress 0.1998022
## Run 69 stress 0.1987399
## ... Procrustes: rmse 0.01450486 max resid 0.1144013
## Run 70 stress 0.1988423
## ... Procrustes: rmse 0.02811877 max resid 0.1712319
## Run 71 stress 0.2019165
## Run 72 stress 0.199738
## Run 73 stress 0.2018132
## Run 74 stress 0.1988518
## ... Procrustes: rmse 0.01858626 max resid 0.1085192
## Run 75 stress 0.1998056
## Run 76 stress 0.2002413
## Run 77 stress 0.19943
## Run 78 stress 0.2025355
## Run 79 stress 0.2032316
## Run 80 stress 0.1988516
## ... Procrustes: rmse 0.03298364 max resid 0.171229
## Run 81 stress 0.2006195
## Run 82 stress 0.2003978
## Run 83 stress 0.1996292
## Run 84 stress 0.2036948
## Run 85 stress 0.2002551
## Run 86 stress 0.1994199
## Run 87 stress 0.2009998
## Run 88 stress 0.199742
## Run 89 stress 0.1989888
## ... Procrustes: rmse 0.02795746 max resid 0.1698642
## Run 90 stress 0.1987502
## ... Procrustes: rmse 0.02321445 max resid 0.1700351
## Run 91 stress 0.1986695
## ... Procrustes: rmse 0.01396068 max resid 0.1183785
## Run 92 stress 0.2031203
## Run 93 stress 0.1996764
## Run 94 stress 0.1986937
## ... Procrustes: rmse 0.02951644 max resid 0.1700355
## Run 95 stress 0.2049831
## Run 96 stress 0.2019955
## Run 97 stress 0.1990337
## Run 98 stress 0.1991023
## Run 99 stress 0.2016464
## Run 100 stress 0.2021279
## Run 101 stress 0.2019234
## Run 102 stress 0.2054546
## Run 103 stress 0.1994465
## Run 104 stress 0.2013532
## Run 105 stress 0.2021416
## Run 106 stress 0.1988531
## ... Procrustes: rmse 0.01992621 max resid 0.1327548
## Run 107 stress 0.2017347
## Run 108 stress 0.2012012
## Run 109 stress 0.1992088
## Run 110 stress 0.2011732
## Run 111 stress 0.2013295
## Run 112 stress 0.1996859
## Run 113 stress 0.1986455
## ... Procrustes: rmse 0.02258753 max resid 0.1713427
## Run 114 stress 0.1989197
## ... Procrustes: rmse 0.01711939 max resid 0.1377771
## Run 115 stress 0.2030923
## Run 116 stress 0.2006961
## Run 117 stress 0.1993259
## Run 118 stress 0.1991269
## Run 119 stress 0.1997224
## Run 120 stress 0.2012424
## Run 121 stress 0.2037182
## Run 122 stress 0.2005638
## Run 123 stress 0.1997293
## Run 124 stress 0.2019241
## Run 125 stress 0.1990247
## Run 126 stress 0.2019549
## Run 127 stress 0.2009573
## Run 128 stress 0.1991389
## Run 129 stress 0.199177
## Run 130 stress 0.2006781
## Run 131 stress 0.1993308
## Run 132 stress 0.2011143
## Run 133 stress 0.1989094
## ... Procrustes: rmse 0.01816284 max resid 0.1119851
## Run 134 stress 0.1990183
## ... Procrustes: rmse 0.02386426 max resid 0.1840146
## Run 135 stress 0.2005903
## Run 136 stress 0.2021391
## Run 137 stress 0.2007415
## Run 138 stress 0.2008611
## Run 139 stress 0.2041761
## Run 140 stress 0.2005297
## Run 141 stress 0.1988056
## ... Procrustes: rmse 0.02600904 max resid 0.17333
## Run 142 stress 0.1991846
## Run 143 stress 0.2029315
## Run 144 stress 0.1992152
## Run 145 stress 0.1987048
## ... Procrustes: rmse 0.02909313 max resid 0.1718795
## Run 146 stress 0.2022298
## Run 147 stress 0.2019817
## Run 148 stress 0.1991584
## Run 149 stress 0.2056617
## Run 150 stress 0.2062946
## Run 151 stress 0.1993116
## Run 152 stress 0.2087897
## Run 153 stress 0.2020314
## Run 154 stress 0.1995646
## Run 155 stress 0.2016056
## Run 156 stress 0.1988324
## ... Procrustes: rmse 0.03236438 max resid 0.1709309
## Run 157 stress 0.2007981
## Run 158 stress 0.1994926
## Run 159 stress 0.2023792
## Run 160 stress 0.2023907
## Run 161 stress 0.2003789
## Run 162 stress 0.1991676
## Run 163 stress 0.1991109
## Run 164 stress 0.2052811
## Run 165 stress 0.1999441
## Run 166 stress 0.1990132
## ... Procrustes: rmse 0.02086043 max resid 0.1299902
## Run 167 stress 0.1993029
## Run 168 stress 0.1987453
## ... Procrustes: rmse 0.01616994 max resid 0.1126899
## Run 169 stress 0.1991857
## Run 170 stress 0.1991895
## Run 171 stress 0.1986268
## ... Procrustes: rmse 0.02376548 max resid 0.1696315
## Run 172 stress 0.2033447
## Run 173 stress 0.2000561
## Run 174 stress 0.1998826
## Run 175 stress 0.2021736
## Run 176 stress 0.2048396
## Run 177 stress 0.2018521
## Run 178 stress 0.1989678
## ... Procrustes: rmse 0.02829059 max resid 0.1704619
## Run 179 stress 0.2006452
## Run 180 stress 0.2013944
## Run 181 stress 0.1992953
## Run 182 stress 0.2007075
## Run 183 stress 0.2051405
## Run 184 stress 0.2072591
## Run 185 stress 0.2014528
## Run 186 stress 0.1991768
## Run 187 stress 0.200469
## Run 188 stress 0.2069287
## Run 189 stress 0.2007361
## Run 190 stress 0.1991609
## Run 191 stress 0.2003093
## Run 192 stress 0.1987035
## ... Procrustes: rmse 0.02781393 max resid 0.1711951
## Run 193 stress 0.1988301
## ... Procrustes: rmse 0.03244099 max resid 0.1706084
## Run 194 stress 0.199728
## Run 195 stress 0.1992211
## Run 196 stress 0.2016071
## Run 197 stress 0.2022122
## Run 198 stress 0.200091
## Run 199 stress 0.1985698
## ... Procrustes: rmse 0.01845712 max resid 0.1712164
## Run 200 stress 0.1992101
## Run 201 stress 0.1997964
## Run 202 stress 0.2044335
## Run 203 stress 0.1988109
## ... Procrustes: rmse 0.02262543 max resid 0.1708955
## Run 204 stress 0.2026284
## Run 205 stress 0.1993863
## Run 206 stress 0.1993545
## Run 207 stress 0.2032916
## Run 208 stress 0.1986746
## ... Procrustes: rmse 0.02752487 max resid 0.1699195
## Run 209 stress 0.2021752
## Run 210 stress 0.1994707
## Run 211 stress 0.2038882
## Run 212 stress 0.1989521
## ... Procrustes: rmse 0.03190345 max resid 0.171515
## Run 213 stress 0.2025844
## Run 214 stress 0.2016695
## Run 215 stress 0.1992049
## Run 216 stress 0.2013576
## Run 217 stress 0.2017951
## Run 218 stress 0.2024938
## Run 219 stress 0.1995634
## Run 220 stress 0.2008625
## Run 221 stress 0.1996773
## Run 222 stress 0.198711
## ... Procrustes: rmse 0.01408558 max resid 0.1189546
## Run 223 stress 0.1984912
## ... New best solution
## ... Procrustes: rmse 0.01773033 max resid 0.1713718
## Run 224 stress 0.2059708
## Run 225 stress 0.2007443
## Run 226 stress 0.2014865
## Run 227 stress 0.2004626
## Run 228 stress 0.199306
## Run 229 stress 0.2023975
## Run 230 stress 0.1985924
## ... Procrustes: rmse 0.01913317 max resid 0.171534
## Run 231 stress 0.2026194
## Run 232 stress 0.1989433
## ... Procrustes: rmse 0.02504959 max resid 0.1685237
## Run 233 stress 0.1996404
## Run 234 stress 0.1987875
## ... Procrustes: rmse 0.0206487 max resid 0.1792166
## Run 235 stress 0.1990171
## Run 236 stress 0.2023326
## Run 237 stress 0.1995038
## Run 238 stress 0.2029404
## Run 239 stress 0.2021112
## Run 240 stress 0.2065854
## Run 241 stress 0.1992441
## Run 242 stress 0.2051455
## Run 243 stress 0.2010685
## Run 244 stress 0.2008611
## Run 245 stress 0.1990731
## Run 246 stress 0.2018966
## Run 247 stress 0.200944
## Run 248 stress 0.2037318
## Run 249 stress 0.2028264
## Run 250 stress 0.2039525
## Run 251 stress 0.2009541
## Run 252 stress 0.2022857
## Run 253 stress 0.1990767
## Run 254 stress 0.2009144
## Run 255 stress 0.1997707
## Run 256 stress 0.199808
## Run 257 stress 0.2026422
## Run 258 stress 0.2014637
## Run 259 stress 0.199019
## Run 260 stress 0.2007184
## Run 261 stress 0.1996912
## Run 262 stress 0.2022896
## Run 263 stress 0.1993412
## Run 264 stress 0.1991661
## Run 265 stress 0.199218
## Run 266 stress 0.2006018
## Run 267 stress 0.2052953
## Run 268 stress 0.2021999
## Run 269 stress 0.2007348
## Run 270 stress 0.1986738
## ... Procrustes: rmse 0.01667522 max resid 0.1672448
## Run 271 stress 0.2028079
## Run 272 stress 0.1992426
## Run 273 stress 0.2071782
## Run 274 stress 0.202121
## Run 275 stress 0.1991092
## Run 276 stress 0.2009635
## Run 277 stress 0.1985208
## ... Procrustes: rmse 0.01761068 max resid 0.171406
## Run 278 stress 0.2006164
## Run 279 stress 0.2038415
## Run 280 stress 0.2030213
## Run 281 stress 0.1994425
## Run 282 stress 0.2023736
## Run 283 stress 0.1994009
## Run 284 stress 0.2025765
## Run 285 stress 0.1991607
## Run 286 stress 0.2042163
## Run 287 stress 0.1996912
## Run 288 stress 0.2026504
## Run 289 stress 0.1991734
## Run 290 stress 0.2014677
## Run 291 stress 0.2014547
## Run 292 stress 0.2037769
## Run 293 stress 0.1991262
## Run 294 stress 0.2006943
## Run 295 stress 0.2017624
## Run 296 stress 0.2068716
## Run 297 stress 0.1986592
## ... Procrustes: rmse 0.01645219 max resid 0.1661095
## Run 298 stress 0.1991897
## Run 299 stress 0.2025279
## Run 300 stress 0.2004124
## Run 301 stress 0.2025142
## Run 302 stress 0.1991954
## Run 303 stress 0.2013155
## Run 304 stress 0.1991675
## Run 305 stress 0.1986413
## ... Procrustes: rmse 0.02030378 max resid 0.1678261
## Run 306 stress 0.1999362
## Run 307 stress 0.2019689
## Run 308 stress 0.2012059
## Run 309 stress 0.2039098
## Run 310 stress 0.1996977
## Run 311 stress 0.2011775
## Run 312 stress 0.1985632
## ... Procrustes: rmse 0.01981513 max resid 0.1721762
## Run 313 stress 0.199368
## Run 314 stress 0.1990472
## Run 315 stress 0.2061646
## Run 316 stress 0.2032448
## Run 317 stress 0.1992963
## Run 318 stress 0.1989293
## ... Procrustes: rmse 0.0271283 max resid 0.1473412
## Run 319 stress 0.208852
## Run 320 stress 0.2003459
## Run 321 stress 0.2080626
## Run 322 stress 0.2017155
## Run 323 stress 0.2007567
## Run 324 stress 0.1988029
## ... Procrustes: rmse 0.01917766 max resid 0.1446588
## Run 325 stress 0.1990415
## Run 326 stress 0.1993851
## Run 327 stress 0.204022
## Run 328 stress 0.1989178
## ... Procrustes: rmse 0.01969818 max resid 0.1675948
## Run 329 stress 0.2026988
## Run 330 stress 0.1991149
## Run 331 stress 0.2032537
## Run 332 stress 0.1991057
## Run 333 stress 0.2044222
## Run 334 stress 0.1992679
## Run 335 stress 0.2032766
## Run 336 stress 0.2003392
## Run 337 stress 0.2013977
## Run 338 stress 0.1989409
## ... Procrustes: rmse 0.02564664 max resid 0.1406273
## Run 339 stress 0.1987608
## ... Procrustes: rmse 0.01456123 max resid 0.1320128
## Run 340 stress 0.1990519
## Run 341 stress 0.2006603
## Run 342 stress 0.203316
## Run 343 stress 0.1991752
## Run 344 stress 0.1999739
## Run 345 stress 0.1988269
## ... Procrustes: rmse 0.0191455 max resid 0.1702959
## Run 346 stress 0.2060452
## Run 347 stress 0.2025583
## Run 348 stress 0.2036732
## Run 349 stress 0.199137
## Run 350 stress 0.2024447
## Run 351 stress 0.2012564
## Run 352 stress 0.1992833
## Run 353 stress 0.2016963
## Run 354 stress 0.1988612
## ... Procrustes: rmse 0.02104022 max resid 0.136064
## Run 355 stress 0.2003867
## Run 356 stress 0.2017984
## Run 357 stress 0.2026786
## Run 358 stress 0.1989599
## ... Procrustes: rmse 0.02595636 max resid 0.1705967
## Run 359 stress 0.2016505
## Run 360 stress 0.2073217
## Run 361 stress 0.1987311
## ... Procrustes: rmse 0.02311892 max resid 0.1677914
## Run 362 stress 0.2021319
## Run 363 stress 0.2001224
## Run 364 stress 0.2037962
## Run 365 stress 0.1992321
## Run 366 stress 0.1988436
## ... Procrustes: rmse 0.02614391 max resid 0.1440492
## Run 367 stress 0.2007991
## Run 368 stress 0.2085751
## Run 369 stress 0.2019143
## Run 370 stress 0.1985995
## ... Procrustes: rmse 0.01861028 max resid 0.1700012
## Run 371 stress 0.1996423
## Run 372 stress 0.2018191
## Run 373 stress 0.200751
## Run 374 stress 0.2062477
## Run 375 stress 0.1986685
## ... Procrustes: rmse 0.01856171 max resid 0.1642752
## Run 376 stress 0.1990513
## Run 377 stress 0.201043
## Run 378 stress 0.1987349
## ... Procrustes: rmse 0.01336608 max resid 0.0959004
## Run 379 stress 0.2055743
## Run 380 stress 0.2022338
## Run 381 stress 0.2078233
## Run 382 stress 0.2003906
## Run 383 stress 0.2021645
## Run 384 stress 0.2053309
## Run 385 stress 0.2022185
## Run 386 stress 0.2004115
## Run 387 stress 0.2010409
## Run 388 stress 0.2052304
## Run 389 stress 0.2027439
## Run 390 stress 0.1994091
## Run 391 stress 0.2048968
## Run 392 stress 0.202398
## Run 393 stress 0.2002044
## Run 394 stress 0.2029158
## Run 395 stress 0.1988468
## ... Procrustes: rmse 0.02602254 max resid 0.1405422
## Run 396 stress 0.19928
## Run 397 stress 0.202913
## Run 398 stress 0.19917
## Run 399 stress 0.1986618
## ... Procrustes: rmse 0.01350005 max resid 0.1330669
## Run 400 stress 0.1990505
## Run 401 stress 0.1997257
## Run 402 stress 0.19916
## Run 403 stress 0.200581
## Run 404 stress 0.2003484
## Run 405 stress 0.2049949
## Run 406 stress 0.1991664
## Run 407 stress 0.2022686
## Run 408 stress 0.201549
## Run 409 stress 0.1986109
## ... Procrustes: rmse 0.01121253 max resid 0.09101894
## Run 410 stress 0.1994125
## Run 411 stress 0.2049857
## Run 412 stress 0.2002931
## Run 413 stress 0.2038996
## Run 414 stress 0.2063077
## Run 415 stress 0.2030575
## Run 416 stress 0.1989817
## ... Procrustes: rmse 0.02316748 max resid 0.1639358
## Run 417 stress 0.2029164
## Run 418 stress 0.1987988
## ... Procrustes: rmse 0.0182961 max resid 0.136394
## Run 419 stress 0.1998191
## Run 420 stress 0.2012493
## Run 421 stress 0.1991032
## Run 422 stress 0.202147
## Run 423 stress 0.1987808
## ... Procrustes: rmse 0.0178711 max resid 0.1380437
## Run 424 stress 0.2023408
## Run 425 stress 0.2058983
## Run 426 stress 0.1991732
## Run 427 stress 0.1989517
## ... Procrustes: rmse 0.01886185 max resid 0.1390814
## Run 428 stress 0.1988258
## ... Procrustes: rmse 0.02129168 max resid 0.1772656
## Run 429 stress 0.2022351
## Run 430 stress 0.2007435
## Run 431 stress 0.20088
## Run 432 stress 0.2004221
## Run 433 stress 0.1987353
## ... Procrustes: rmse 0.02384011 max resid 0.1438022
## Run 434 stress 0.1996529
## Run 435 stress 0.2026357
## Run 436 stress 0.2043552
## Run 437 stress 0.2023713
## Run 438 stress 0.2003803
## Run 439 stress 0.1989239
## ... Procrustes: rmse 0.02680498 max resid 0.1425762
## Run 440 stress 0.2002909
## Run 441 stress 0.2007285
## Run 442 stress 0.2019365
## Run 443 stress 0.2001735
## Run 444 stress 0.2003657
## Run 445 stress 0.1986067
## ... Procrustes: rmse 0.007991469 max resid 0.04801641
## Run 446 stress 0.1989321
## ... Procrustes: rmse 0.0159085 max resid 0.1351729
## Run 447 stress 0.2057908
## Run 448 stress 0.1986303
## ... Procrustes: rmse 0.0151653 max resid 0.1387526
## Run 449 stress 0.1995467
## Run 450 stress 0.2065196
## Run 451 stress 0.1986612
## ... Procrustes: rmse 0.0223294 max resid 0.1426928
## Run 452 stress 0.2004796
## Run 453 stress 0.1995172
## Run 454 stress 0.1988918
## ... Procrustes: rmse 0.02700406 max resid 0.1433678
## Run 455 stress 0.2005308
## Run 456 stress 0.1990586
## Run 457 stress 0.1992004
## Run 458 stress 0.1996166
## Run 459 stress 0.1998616
## Run 460 stress 0.2040413
## Run 461 stress 0.2042319
## Run 462 stress 0.2019378
## Run 463 stress 0.1986483
## ... Procrustes: rmse 0.01184867 max resid 0.09592897
## Run 464 stress 0.1986762
## ... Procrustes: rmse 0.0217787 max resid 0.1417871
## Run 465 stress 0.1997891
## Run 466 stress 0.1993738
## Run 467 stress 0.201443
## Run 468 stress 0.2027473
## Run 469 stress 0.1992299
## Run 470 stress 0.1998202
## Run 471 stress 0.1991743
## Run 472 stress 0.1991663
## Run 473 stress 0.2015746
## Run 474 stress 0.1994073
## Run 475 stress 0.1993344
## Run 476 stress 0.2014867
## Run 477 stress 0.2066334
## Run 478 stress 0.1987782
## ... Procrustes: rmse 0.01753591 max resid 0.1355994
## Run 479 stress 0.200297
## Run 480 stress 0.2021401
## Run 481 stress 0.2048378
## Run 482 stress 0.2065698
## Run 483 stress 0.2007548
## Run 484 stress 0.2031585
## Run 485 stress 0.2100238
## Run 486 stress 0.199304
## Run 487 stress 0.1993894
## Run 488 stress 0.1991846
## Run 489 stress 0.2011221
## Run 490 stress 0.2022824
## Run 491 stress 0.2033317
## Run 492 stress 0.1991575
## Run 493 stress 0.2087248
## Run 494 stress 0.1988533
## ... Procrustes: rmse 0.02039052 max resid 0.1379401
## Run 495 stress 0.2046087
## Run 496 stress 0.1986408
## ... Procrustes: rmse 0.02043356 max resid 0.1687092
## Run 497 stress 0.1992631
## Run 498 stress 0.2015601
## Run 499 stress 0.2019134
## Run 500 stress 0.2054894
## Run 501 stress 0.2002991
## Run 502 stress 0.1993832
## Run 503 stress 0.2000892
## Run 504 stress 0.2039328
## Run 505 stress 0.1990072
## Run 506 stress 0.1986732
## ... Procrustes: rmse 0.0104371 max resid 0.1014971
## Run 507 stress 0.1991884
## Run 508 stress 0.2017655
## Run 509 stress 0.1989721
## ... Procrustes: rmse 0.02338846 max resid 0.165625
## Run 510 stress 0.2001217
## Run 511 stress 0.1991614
## Run 512 stress 0.2027388
## Run 513 stress 0.1991206
## Run 514 stress 0.2041714
## Run 515 stress 0.2047861
## Run 516 stress 0.205032
## Run 517 stress 0.2007435
## Run 518 stress 0.2003388
## Run 519 stress 0.1990681
## Run 520 stress 0.1988504
## ... Procrustes: rmse 0.02125788 max resid 0.1796563
## Run 521 stress 0.1991786
## Run 522 stress 0.199755
## Run 523 stress 0.2011018
## Run 524 stress 0.1986369
## ... Procrustes: rmse 0.01964026 max resid 0.1689584
## Run 525 stress 0.2009477
## Run 526 stress 0.1992634
## Run 527 stress 0.2020909
## Run 528 stress 0.2011073
## Run 529 stress 0.200599
## Run 530 stress 0.1989779
## ... Procrustes: rmse 0.02253931 max resid 0.1388132
## Run 531 stress 0.2069369
## Run 532 stress 0.1998343
## Run 533 stress 0.1991122
## Run 534 stress 0.2018708
## Run 535 stress 0.1986035
## ... Procrustes: rmse 0.01854648 max resid 0.1691484
## Run 536 stress 0.1991606
## Run 537 stress 0.2028867
## Run 538 stress 0.1985637
## ... Procrustes: rmse 0.00495467 max resid 0.03479443
## Run 539 stress 0.2001315
## Run 540 stress 0.1995799
## Run 541 stress 0.2076461
## Run 542 stress 0.2065378
## Run 543 stress 0.1989535
## ... Procrustes: rmse 0.02524189 max resid 0.1415103
## Run 544 stress 0.2018671
## Run 545 stress 0.1990145
## Run 546 stress 0.1987056
## ... Procrustes: rmse 0.02118785 max resid 0.1472073
## Run 547 stress 0.2028308
## Run 548 stress 0.1992645
## Run 549 stress 0.2025836
## Run 550 stress 0.1991243
## Run 551 stress 0.1986862
## ... Procrustes: rmse 0.01647469 max resid 0.150392
## Run 552 stress 0.2044257
## Run 553 stress 0.2041052
## Run 554 stress 0.2024971
## Run 555 stress 0.2025812
## Run 556 stress 0.2038905
## Run 557 stress 0.1991462
## Run 558 stress 0.1997909
## Run 559 stress 0.2036684
## Run 560 stress 0.1988405
## ... Procrustes: rmse 0.02414166 max resid 0.1406465
## Run 561 stress 0.2049576
## Run 562 stress 0.1991416
## Run 563 stress 0.2020148
## Run 564 stress 0.2047058
## Run 565 stress 0.2015827
## Run 566 stress 0.2008732
## Run 567 stress 0.2013834
## Run 568 stress 0.2024219
## Run 569 stress 0.1992404
## Run 570 stress 0.2019534
## Run 571 stress 0.199672
## Run 572 stress 0.20154
## Run 573 stress 0.1996973
## Run 574 stress 0.1986907
## ... Procrustes: rmse 0.02156798 max resid 0.1681474
## Run 575 stress 0.2000584
## Run 576 stress 0.1991
## Run 577 stress 0.2022946
## Run 578 stress 0.2002615
## Run 579 stress 0.2027747
## Run 580 stress 0.1994205
## Run 581 stress 0.1986019
## ... Procrustes: rmse 0.01110782 max resid 0.09322401
## Run 582 stress 0.1997304
## Run 583 stress 0.2028305
## Run 584 stress 0.2015479
## Run 585 stress 0.1986946
## ... Procrustes: rmse 0.02047597 max resid 0.139702
## Run 586 stress 0.2034813
## Run 587 stress 0.2004891
## Run 588 stress 0.1994006
## Run 589 stress 0.2016301
## Run 590 stress 0.1989248
## ... Procrustes: rmse 0.02676526 max resid 0.1430197
## Run 591 stress 0.2007402
## Run 592 stress 0.1987863
## ... Procrustes: rmse 0.01997184 max resid 0.1622066
## Run 593 stress 0.2036853
## Run 594 stress 0.2034341
## Run 595 stress 0.1997136
## Run 596 stress 0.2048613
## Run 597 stress 0.2024261
## Run 598 stress 0.1992226
## Run 599 stress 0.1988838
## ... Procrustes: rmse 0.01710005 max resid 0.135213
## Run 600 stress 0.1988602
## ... Procrustes: rmse 0.01880647 max resid 0.171163
## Run 601 stress 0.203591
## Run 602 stress 0.1987943
## ... Procrustes: rmse 0.02168424 max resid 0.1437009
## Run 603 stress 0.2025379
## Run 604 stress 0.2034697
## Run 605 stress 0.201036
## Run 606 stress 0.2016669
## Run 607 stress 0.2021754
## Run 608 stress 0.2044511
## Run 609 stress 0.2006244
## Run 610 stress 0.1997862
## Run 611 stress 0.199041
## Run 612 stress 0.1993785
## Run 613 stress 0.19974
## Run 614 stress 0.1989058
## ... Procrustes: rmse 0.02310342 max resid 0.1454902
## Run 615 stress 0.1988474
## ... Procrustes: rmse 0.02259532 max resid 0.1416227
## Run 616 stress 0.198858
## ... Procrustes: rmse 0.0228228 max resid 0.1625409
## Run 617 stress 0.2037488
## Run 618 stress 0.2013586
## Run 619 stress 0.2042178
## Run 620 stress 0.2001878
## Run 621 stress 0.2026393
## Run 622 stress 0.2004966
## Run 623 stress 0.1991642
## Run 624 stress 0.2039664
## Run 625 stress 0.201141
## Run 626 stress 0.1991785
## Run 627 stress 0.1998065
## Run 628 stress 0.2000673
## Run 629 stress 0.1993308
## Run 630 stress 0.2013287
## Run 631 stress 0.1988618
## ... Procrustes: rmse 0.02497703 max resid 0.1387103
## Run 632 stress 0.204586
## Run 633 stress 0.2020663
## Run 634 stress 0.1996984
## Run 635 stress 0.1993298
## Run 636 stress 0.2024929
## Run 637 stress 0.1989684
## ... Procrustes: rmse 0.02756382 max resid 0.1454735
## Run 638 stress 0.2020633
## Run 639 stress 0.2045711
## Run 640 stress 0.2008809
## Run 641 stress 0.1990141
## Run 642 stress 0.1990122
## Run 643 stress 0.1997368
## Run 644 stress 0.1986738
## ... Procrustes: rmse 0.0194998 max resid 0.1663586
## Run 645 stress 0.1985882
## ... Procrustes: rmse 0.01838105 max resid 0.1671357
## Run 646 stress 0.2003448
## Run 647 stress 0.2061547
## Run 648 stress 0.1988823
## ... Procrustes: rmse 0.02684273 max resid 0.142012
## Run 649 stress 0.2023909
## Run 650 stress 0.2018409
## Run 651 stress 0.1988065
## ... Procrustes: rmse 0.01830718 max resid 0.1462947
## Run 652 stress 0.1988909
## ... Procrustes: rmse 0.01694858 max resid 0.1366914
## Run 653 stress 0.1993351
## Run 654 stress 0.2005719
## Run 655 stress 0.2013265
## Run 656 stress 0.2039915
## Run 657 stress 0.2072421
## Run 658 stress 0.1997402
## Run 659 stress 0.204453
## Run 660 stress 0.1987989
## ... Procrustes: rmse 0.01822103 max resid 0.135578
## Run 661 stress 0.2053136
## Run 662 stress 0.1991592
## Run 663 stress 0.198875
## ... Procrustes: rmse 0.02095661 max resid 0.1414295
## Run 664 stress 0.2002318
## Run 665 stress 0.2008881
## Run 666 stress 0.2000439
## Run 667 stress 0.1988102
## ... Procrustes: rmse 0.01849812 max resid 0.1373347
## Run 668 stress 0.1987354
## ... Procrustes: rmse 0.01091993 max resid 0.06499744
## Run 669 stress 0.1994997
## Run 670 stress 0.2029557
## Run 671 stress 0.2024706
## Run 672 stress 0.2026116
## Run 673 stress 0.1992405
## Run 674 stress 0.2014647
## Run 675 stress 0.1993338
## Run 676 stress 0.2020824
## Run 677 stress 0.19867
## ... Procrustes: rmse 0.02274431 max resid 0.1448562
## Run 678 stress 0.2048384
## Run 679 stress 0.200284
## Run 680 stress 0.2002257
## Run 681 stress 0.1989887
## ... Procrustes: rmse 0.02677931 max resid 0.1677287
## Run 682 stress 0.2083509
## Run 683 stress 0.1986493
## ... Procrustes: rmse 0.008871632 max resid 0.04943209
## Run 684 stress 0.1988307
## ... Procrustes: rmse 0.01449009 max resid 0.1015514
## Run 685 stress 0.2020483
## Run 686 stress 0.2022033
## Run 687 stress 0.2053077
## Run 688 stress 0.2026392
## Run 689 stress 0.1998337
## Run 690 stress 0.1997464
## Run 691 stress 0.1990484
## Run 692 stress 0.199021
## Run 693 stress 0.2012827
## Run 694 stress 0.2041858
## Run 695 stress 0.2018441
## Run 696 stress 0.2045235
## Run 697 stress 0.2021418
## Run 698 stress 0.2020666
## Run 699 stress 0.2044883
## Run 700 stress 0.1991144
## Run 701 stress 0.2017284
## Run 702 stress 0.19918
## Run 703 stress 0.2004179
## Run 704 stress 0.202071
## Run 705 stress 0.1987644
## ... Procrustes: rmse 0.02050274 max resid 0.1825274
## Run 706 stress 0.1990239
## Run 707 stress 0.2073208
## Run 708 stress 0.1990182
## Run 709 stress 0.1986668
## ... Procrustes: rmse 0.01933432 max resid 0.1665206
## Run 710 stress 0.2027118
## Run 711 stress 0.2038226
## Run 712 stress 0.2058146
## Run 713 stress 0.2005921
## Run 714 stress 0.2022129
## Run 715 stress 0.2026442
## Run 716 stress 0.2008159
## Run 717 stress 0.2037802
## Run 718 stress 0.2019701
## Run 719 stress 0.1998588
## Run 720 stress 0.2028501
## Run 721 stress 0.2023026
## Run 722 stress 0.2064202
## Run 723 stress 0.199065
## Run 724 stress 0.1991403
## Run 725 stress 0.2015405
## Run 726 stress 0.2007991
## Run 727 stress 0.1989076
## ... Procrustes: rmse 0.02512583 max resid 0.1818442
## Run 728 stress 0.2021169
## Run 729 stress 0.1991388
## Run 730 stress 0.2004083
## Run 731 stress 0.2046412
## Run 732 stress 0.2051396
## Run 733 stress 0.1993771
## Run 734 stress 0.1992853
## Run 735 stress 0.2020609
## Run 736 stress 0.1991536
## Run 737 stress 0.2050524
## Run 738 stress 0.2032992
## Run 739 stress 0.2012531
## Run 740 stress 0.1988448
## ... Procrustes: rmse 0.02562965 max resid 0.1405705
## Run 741 stress 0.2010867
## Run 742 stress 0.1992258
## Run 743 stress 0.1987261
## ... Procrustes: rmse 0.01485932 max resid 0.1347667
## Run 744 stress 0.2022138
## Run 745 stress 0.2003978
## Run 746 stress 0.1992092
## Run 747 stress 0.2030846
## Run 748 stress 0.2003135
## Run 749 stress 0.2007524
## Run 750 stress 0.2014736
## Run 751 stress 0.2008572
## Run 752 stress 0.2013288
## Run 753 stress 0.2017241
## Run 754 stress 0.2039402
## Run 755 stress 0.1996482
## Run 756 stress 0.2002254
## Run 757 stress 0.2024773
## Run 758 stress 0.1991269
## Run 759 stress 0.2003948
## Run 760 stress 0.2000613
## Run 761 stress 0.1987633
## ... Procrustes: rmse 0.02339543 max resid 0.1662057
## Run 762 stress 0.1985619
## ... Procrustes: rmse 0.01827523 max resid 0.1709286
## Run 763 stress 0.2018599
## Run 764 stress 0.2068496
## Run 765 stress 0.1990354
## Run 766 stress 0.1991407
## Run 767 stress 0.2011272
## Run 768 stress 0.2034131
## Run 769 stress 0.1993313
## Run 770 stress 0.2019801
## Run 771 stress 0.1990974
## Run 772 stress 0.2003637
## Run 773 stress 0.201668
## Run 774 stress 0.2026323
## Run 775 stress 0.1991242
## Run 776 stress 0.2015436
## Run 777 stress 0.2034487
## Run 778 stress 0.2015309
## Run 779 stress 0.202289
## Run 780 stress 0.1986747
## ... Procrustes: rmse 0.02009041 max resid 0.1428834
## Run 781 stress 0.1991187
## Run 782 stress 0.2020734
## Run 783 stress 0.2047688
## Run 784 stress 0.1990142
## Run 785 stress 0.2037572
## Run 786 stress 0.1991312
## Run 787 stress 0.1991214
## Run 788 stress 0.202534
## Run 789 stress 0.202178
## Run 790 stress 0.1992981
## Run 791 stress 0.2025517
## Run 792 stress 0.2036247
## Run 793 stress 0.1986995
## ... Procrustes: rmse 0.0209905 max resid 0.1457293
## Run 794 stress 0.2001074
## Run 795 stress 0.199178
## Run 796 stress 0.2007435
## Run 797 stress 0.2014563
## Run 798 stress 0.2031501
## Run 799 stress 0.19915
## Run 800 stress 0.2005952
## Run 801 stress 0.1988503
## ... Procrustes: rmse 0.01821611 max resid 0.1684881
## Run 802 stress 0.2034334
## Run 803 stress 0.1987637
## ... Procrustes: rmse 0.02138735 max resid 0.1429999
## Run 804 stress 0.2011508
## Run 805 stress 0.2029092
## Run 806 stress 0.1986214
## ... Procrustes: rmse 0.01893333 max resid 0.169673
## Run 807 stress 0.2050758
## Run 808 stress 0.1991613
## Run 809 stress 0.2032006
## Run 810 stress 0.2009335
## Run 811 stress 0.2014617
## Run 812 stress 0.2020048
## Run 813 stress 0.2004156
## Run 814 stress 0.1992308
## Run 815 stress 0.1986927
## ... Procrustes: rmse 0.01683715 max resid 0.1474331
## Run 816 stress 0.1986993
## ... Procrustes: rmse 0.01930659 max resid 0.1369754
## Run 817 stress 0.2005487
## Run 818 stress 0.2021351
## Run 819 stress 0.2002859
## Run 820 stress 0.2055323
## Run 821 stress 0.1997226
## Run 822 stress 0.1989517
## ... Procrustes: rmse 0.02052125 max resid 0.1640415
## Run 823 stress 0.2009705
## Run 824 stress 0.2036171
## Run 825 stress 0.2010037
## Run 826 stress 0.2017672
## Run 827 stress 0.2039537
## Run 828 stress 0.2025844
## Run 829 stress 0.2026681
## Run 830 stress 0.1999901
## Run 831 stress 0.2067949
## Run 832 stress 0.1988937
## ... Procrustes: rmse 0.01976112 max resid 0.1331059
## Run 833 stress 0.2067978
## Run 834 stress 0.2025308
## Run 835 stress 0.2009436
## Run 836 stress 0.1990192
## Run 837 stress 0.2006654
## Run 838 stress 0.2043557
## Run 839 stress 0.204113
## Run 840 stress 0.2019374
## Run 841 stress 0.1991222
## Run 842 stress 0.2020171
## Run 843 stress 0.2029479
## Run 844 stress 0.1993644
## Run 845 stress 0.1991044
## Run 846 stress 0.201873
## Run 847 stress 0.2019889
## Run 848 stress 0.2021321
## Run 849 stress 0.2009921
## Run 850 stress 0.1990226
## Run 851 stress 0.205095
## Run 852 stress 0.2000952
## Run 853 stress 0.1990091
## Run 854 stress 0.20203
## Run 855 stress 0.2045591
## Run 856 stress 0.2063274
## Run 857 stress 0.1993386
## Run 858 stress 0.2018171
## Run 859 stress 0.2002256
## Run 860 stress 0.2021338
## Run 861 stress 0.2010103
## Run 862 stress 0.1994052
## Run 863 stress 0.2001932
## Run 864 stress 0.2002461
## Run 865 stress 0.1989542
## ... Procrustes: rmse 0.02632591 max resid 0.1413673
## Run 866 stress 0.2028376
## Run 867 stress 0.2032624
## Run 868 stress 0.2011599
## Run 869 stress 0.2038658
## Run 870 stress 0.2005662
## Run 871 stress 0.1994379
## Run 872 stress 0.1986297
## ... Procrustes: rmse 0.01492073 max resid 0.137669
## Run 873 stress 0.1989731
## ... Procrustes: rmse 0.02008362 max resid 0.1285267
## Run 874 stress 0.1998107
## Run 875 stress 0.2002326
## Run 876 stress 0.1990943
## Run 877 stress 0.1993051
## Run 878 stress 0.1991943
## Run 879 stress 0.2016233
## Run 880 stress 0.2019455
## Run 881 stress 0.1998036
## Run 882 stress 0.2066104
## Run 883 stress 0.2025454
## Run 884 stress 0.1988383
## ... Procrustes: rmse 0.02159349 max resid 0.1829847
## Run 885 stress 0.2014614
## Run 886 stress 0.1998104
## Run 887 stress 0.1987441
## ... Procrustes: rmse 0.01995285 max resid 0.1662434
## Run 888 stress 0.1990352
## Run 889 stress 0.2001747
## Run 890 stress 0.2006182
## Run 891 stress 0.1990192
## Run 892 stress 0.2013548
## Run 893 stress 0.2037468
## Run 894 stress 0.1994997
## Run 895 stress 0.1989713
## ... Procrustes: rmse 0.0217352 max resid 0.1382353
## Run 896 stress 0.2081567
## Run 897 stress 0.1986517
## ... Procrustes: rmse 0.01201379 max resid 0.08771099
## Run 898 stress 0.2034279
## Run 899 stress 0.2030125
## Run 900 stress 0.1991215
## Run 901 stress 0.2000278
## Run 902 stress 0.2045228
## Run 903 stress 0.1994123
## Run 904 stress 0.1988019
## ... Procrustes: rmse 0.02477634 max resid 0.1674513
## Run 905 stress 0.1993287
## Run 906 stress 0.2000729
## Run 907 stress 0.2046898
## Run 908 stress 0.1990183
## Run 909 stress 0.1999495
## Run 910 stress 0.2011285
## Run 911 stress 0.2006487
## Run 912 stress 0.1990014
## Run 913 stress 0.2008053
## Run 914 stress 0.2003953
## Run 915 stress 0.2019991
## Run 916 stress 0.2034195
## Run 917 stress 0.1986585
## ... Procrustes: rmse 0.01262861 max resid 0.09499231
## Run 918 stress 0.2025462
## Run 919 stress 0.2029002
## Run 920 stress 0.2049672
## Run 921 stress 0.2003703
## Run 922 stress 0.2001817
## Run 923 stress 0.2053881
## Run 924 stress 0.2000061
## Run 925 stress 0.2054659
## Run 926 stress 0.1986929
## ... Procrustes: rmse 0.02171009 max resid 0.1684794
## Run 927 stress 0.1996843
## Run 928 stress 0.1991611
## Run 929 stress 0.2016914
## Run 930 stress 0.1992535
## Run 931 stress 0.203878
## Run 932 stress 0.2003175
## Run 933 stress 0.1996764
## Run 934 stress 0.1989184
## ... Procrustes: rmse 0.02590217 max resid 0.1422541
## Run 935 stress 0.1990181
## Run 936 stress 0.1987105
## ... Procrustes: rmse 0.01698843 max resid 0.1672901
## Run 937 stress 0.1991844
## Run 938 stress 0.2009748
## Run 939 stress 0.2019531
## Run 940 stress 0.1997282
## Run 941 stress 0.1991391
## Run 942 stress 0.2007545
## Run 943 stress 0.2009355
## Run 944 stress 0.2033463
## Run 945 stress 0.1989401
## ... Procrustes: rmse 0.02561098 max resid 0.1407032
## Run 946 stress 0.2014832
## Run 947 stress 0.2035425
## Run 948 stress 0.2013074
## Run 949 stress 0.2031958
## Run 950 stress 0.1993708
## Run 951 stress 0.2007878
## Run 952 stress 0.2032039
## Run 953 stress 0.2026133
## Run 954 stress 0.2028267
## Run 955 stress 0.2057689
## Run 956 stress 0.2035232
## Run 957 stress 0.2025934
## Run 958 stress 0.1988262
## ... Procrustes: rmse 0.02530351 max resid 0.1420288
## Run 959 stress 0.1989212
## ... Procrustes: rmse 0.02667251 max resid 0.1426737
## Run 960 stress 0.2015516
## Run 961 stress 0.200508
## Run 962 stress 0.1992802
## Run 963 stress 0.1991276
## Run 964 stress 0.200881
## Run 965 stress 0.2005723
## Run 966 stress 0.1997366
## Run 967 stress 0.2015702
## Run 968 stress 0.2027437
## Run 969 stress 0.2002148
## Run 970 stress 0.1991504
## Run 971 stress 0.1989312
## ... Procrustes: rmse 0.02666753 max resid 0.1428223
## Run 972 stress 0.2031014
## Run 973 stress 0.204282
## Run 974 stress 0.1987996
## ... Procrustes: rmse 0.01163665 max resid 0.0772431
## Run 975 stress 0.1993205
## Run 976 stress 0.1986249
## ... Procrustes: rmse 0.01444663 max resid 0.1370688
## Run 977 stress 0.1995497
## Run 978 stress 0.1995783
## Run 979 stress 0.2020041
## Run 980 stress 0.2036218
## Run 981 stress 0.2035619
## Run 982 stress 0.1991089
## Run 983 stress 0.2056828
## Run 984 stress 0.2051017
## Run 985 stress 0.2054828
## Run 986 stress 0.2032224
## Run 987 stress 0.1992158
## Run 988 stress 0.1989906
## ... Procrustes: rmse 0.02656981 max resid 0.1787974
## Run 989 stress 0.2040098
## Run 990 stress 0.1992773
## Run 991 stress 0.2030067
## Run 992 stress 0.2023959
## Run 993 stress 0.1988437
## ... Procrustes: rmse 0.01620768 max resid 0.1206213
## Run 994 stress 0.2010593
## Run 995 stress 0.1991431
## Run 996 stress 0.2071067
## Run 997 stress 0.1993893
## Run 998 stress 0.204968
## Run 999 stress 0.2053016
## Run 1000 stress 0.2009163
## *** No convergence -- monoMDS stopping criteria:
## 204: no. of iterations >= maxit
## 796: stress ratio > sratmax
#prepare NMDS points and values for graphing
md.nmds.points <- md.nmds$points %>% # take the NMDS points
data.frame() %>% # make a data frame
rownames_to_column(var = "mysamples") %>%
mutate(SampleID = mysamples) %>%
column_to_rownames(var = "mysamples")
meta <- input_toad$map_loaded %>%
rownames_to_column("SampleID")
md.nmds.metadata <- inner_join(x = md.nmds.points, y = meta,
by = "SampleID") %>%
group_by(Life_Stage_Simplified) %>%
ungroup() %>%
dplyr::select(SampleID, MDS1, MDS2, everything()) # select these columns only
md.nmds.metadata.unq <- md.nmds.metadata %>%
dplyr::select(Life_Stage_Simplified) %>%
unique()
# get stress value
stress.md = paste("stress =", round(md.nmds$stress, digits = 4))
# get NMDS hulls
group.chulls <- plyr::ddply(md.nmds.metadata, "Life_Stage_Simplified",
function(df) df[chull(df$MDS1, df$MDS2), ])
# make color pallete
samptype_cols_Toad <- c("#66A61E",
"#E7298A", "#7570B3",
"#D95F02", "#1B9E77")
# Plot NMDS
nmds.plot.toad <- ggplot(md.nmds.metadata, aes(x = MDS1, y = MDS2,
color = Life_Stage_Simplified)) +
geom_polygon(data = group.chulls,
aes(fill = Life_Stage_Simplified),
alpha = 0.15, linetype = 0) +
geom_point(size = 3) +
annotate("text", x = Inf, y = Inf, label = stress.md, hjust = 1, vjust = 1) +
ggtitle("D") +
labs(colour = "Life Stage") +
scale_colour_manual(values = samptype_cols_Toad) +
scale_fill_manual(values = samptype_cols_Toad) +
guides(fill = F, color = guide_legend(order = 1,
override.aes = list(shape = 15,
size = 6))) +
theme(axis.title = element_text(size = 16),
axis.text = element_text(size = 14),
legend.title = element_text(size = 16),
legend.text = element_text(size = 14),
plot.title = element_text(size = 30))
# Make dispersion graph
# PERMDisp on sample type, strata=site
#create B-C species abundance matrix
input_all_rel <- convert_to_relative_abundances(input_all)
BC <- calc_dm(input_all_rel$data_loaded, method = "bray")
life.disper.lump <- betadisper(BC, input_all$map_loaded$Life_Stage_Simplified)
permutest(life.disper.lump)
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 6 0.26325 0.043875 8.5606 999 0.001 ***
## Residuals 162 0.83029 0.005125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
PermDisp_life_lump <- TukeyHSD(life.disper.lump) #Tukey post-hoc
kable(PermDisp_life_lump$group) %>%
kable_styling(bootstrap_options = c("striped","bordered"),
full_width = F, position = "left")| diff | lwr | upr | p adj | |
|---|---|---|---|---|
| Water-Sediment | 0.0524077 | -0.0117179 | 0.1165333 | 0.1889160 |
| Eggs-Sediment | 0.0000077 | -0.0750595 | 0.0750750 | 1.0000000 |
| Tadpole-Sediment | 0.0336625 | -0.0165358 | 0.0838608 | 0.4173172 |
| Metamorph-Sediment | -0.0113815 | -0.0718396 | 0.0490767 | 0.9977330 |
| Subadult-Sediment | -0.0917708 | -0.1964875 | 0.0129459 | 0.1278208 |
| Adult-Sediment | -0.0781873 | -0.1466213 | -0.0097533 | 0.0140724 |
| Eggs-Water | -0.0524000 | -0.1304512 | 0.0256512 | 0.4158847 |
| Tadpole-Water | -0.0187452 | -0.0733049 | 0.0358146 | 0.9473845 |
| Metamorph-Water | -0.0637891 | -0.1279147 | 0.0003365 | 0.0521973 |
| Subadult-Water | -0.1441785 | -0.2510545 | -0.0373025 | 0.0016470 |
| Adult-Water | -0.1305950 | -0.2022896 | -0.0589004 | 0.0000040 |
| Tadpole-Eggs | 0.0336548 | -0.0334255 | 0.1007351 | 0.7457009 |
| Metamorph-Eggs | -0.0113892 | -0.0864564 | 0.0636780 | 0.9993322 |
| Subadult-Eggs | -0.0917785 | -0.2055567 | 0.0219997 | 0.2018068 |
| Adult-Eggs | -0.0781950 | -0.1598229 | 0.0034329 | 0.0700374 |
| Metamorph-Tadpole | -0.0450440 | -0.0952423 | 0.0051543 | 0.1102650 |
| Subadult-Tadpole | -0.1254333 | -0.2245810 | -0.0262857 | 0.0041018 |
| Adult-Tadpole | -0.1118498 | -0.1714140 | -0.0522856 | 0.0000018 |
| Subadult-Metamorph | -0.0803893 | -0.1851060 | 0.0243273 | 0.2542094 |
| Adult-Metamorph | -0.0668058 | -0.1352399 | 0.0016282 | 0.0606342 |
| Adult-Subadult | 0.0135835 | -0.0959318 | 0.1230988 | 0.9997915 |
# write.table(PermDisp_life_lump$group,
# paste0(work_dir,"Figs_Tables/PermDisp_life_lump.txt"),
# sep = "\t")
Disper1 <- data.frame(life.disper.lump$distances)
colnames(Disper1) <- "dispers"
meta <- meta <- input_all$map_loaded %>%
rownames_to_column("SampleID")
Disper2 <- Disper1 %>%
rownames_to_column("SampleID") %>%
inner_join(y = meta,
by = "SampleID")
disper_bars <- ggplot(data = Disper2, aes(x = Life_Stage_Simplified,
y = dispers)) +
geom_boxplot(aes(fill = Life_Stage_Simplified), alpha = 0.6,
outlier.shape = NA) +
ylab("Dispersion") + xlab("") + ggtitle("B") +
geom_point(aes(fill = Life_Stage_Simplified), size = 3, shape = 21, position = position_jitter(
width = 0.25
)) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, colour = "black"),
axis.text.y = element_text(colour = "black"),
axis.title = element_text(size = 16),
axis.text = element_text(size = 14),
legend.title = element_text(size = 16),
legend.text = element_text(size = 14),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
plot.title = element_text(size = 30)) +
labs(fill = "Sample Type") +
scale_fill_manual(values = samptype_cols)
legend <- arrangeGrob(get_legend(disper_bars))
disper_bars_noleg <- disper_bars + theme(legend.position = "none")
# Put all into paneled graph with annotations and shared legend
diversity_fig <- ggpubr::ggarrange(arrangeGrob(OTUrich, disper_bars_noleg),
legend,
arrangeGrob(nmds.plot, nmds.plot.toad),
ncol = 3, widths = c(5,1.5,6))
diversity_fig
# ggsave(paste0(work_dir,"Figs_Tables/FungDivFig.png"), diversity_fig,
# dpi = 800,
# width = 13, height = 8)Permanova for the above plot
# PERMANOVA on sample type, strata=site
Perm <- adonis(BC ~ input_all$map_loaded$Type,
strata = input_all$map_loaded$Location)
kable(Perm$aov.tab) %>%
kable_styling(bootstrap_options = c("striped","bordered"),
full_width = F, position = "left")| Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(>F) | |
|---|---|---|---|---|---|---|
| input_all\(map_loaded\)Type | 1 | 2.035769 | 2.0357689 | 4.870024 | 0.0283355 | 0.001 |
| Residuals | 167 | 69.809399 | 0.4180204 | NA | 0.9716645 | NA |
| Total | 168 | 71.845168 | NA | NA | 1.0000000 | NA |
# PERMANOVA on life stage, strata-site
input_toad_rel <- convert_to_relative_abundances(input_toad)
BC2 <- calc_dm(input_toad_rel$data_loaded, method = "bray")
Perm2 <- adonis(BC2 ~ input_toad$map_loaded$Life_Stage_Simplified,
strata = input_toad$map_loaded$Location)
kable(Perm2$aov.tab) %>%
kable_styling(bootstrap_options = c("striped","bordered"),
full_width = F, position = "left")| Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(>F) | |
|---|---|---|---|---|---|---|
| input_toad\(map_loaded\)Life_Stage_Simplified | 4 | 4.424554 | 1.1061384 | 2.842999 | 0.0872273 | 0.001 |
| Residuals | 119 | 46.299860 | 0.3890745 | NA | 0.9127727 | NA |
| Total | 123 | 50.724413 | NA | NA | 1.0000000 | NA |
# PERMANOVA on all samples by sample type, strata=site
Perm <- adonis(BC ~ input_all$map_loaded$Life_Stage_Simplified,
strata = input_all$map_loaded$Location)
kable(Perm$aov.tab) %>%
kable_styling(bootstrap_options = c("striped","bordered"),
full_width = F, position = "left")| Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(>F) | |
|---|---|---|---|---|---|---|
| input_all\(map_loaded\)Life_Stage_Simplified | 6 | 7.751078 | 1.2918463 | 3.265186 | 0.1078859 | 0.001 |
| Residuals | 162 | 64.094090 | 0.3956425 | NA | 0.8921141 | NA |
| Total | 168 | 71.845168 | NA | NA | 1.0000000 | NA |
# also show the pairwise PERMANOVA
pw_perm <- read.delim(paste0(work_dir,"Figs_Tables/PW_PERM.txt"))
kable(pw_perm) %>%
kable_styling(bootstrap_options = c("striped","bordered"),
full_width = F, position = "left")| X1 | X2 | R2 | pval | pvalBon | pvalFDR |
|---|---|---|---|---|---|
| Water | Sediment | 0.0676608 | 0.001 | 0.078 | 0.008 |
| Water | Eggs.11.15 | 0.0682710 | 0.001 | 0.078 | 0.008 |
| Water | Metamorph.40.46 | 0.0524184 | 0.001 | 0.078 | 0.008 |
| Water | Subadult | 0.0759918 | 0.001 | 0.078 | 0.008 |
| Water | Adult | 0.0945756 | 0.001 | 0.078 | 0.008 |
| Sediment | Eggs.11.15 | 0.0774906 | 0.001 | 0.078 | 0.008 |
| Sediment | Tadpole.25.27 | 0.0828321 | 0.001 | 0.078 | 0.008 |
| Sediment | Tadpole.27.29 | 0.0722623 | 0.001 | 0.078 | 0.008 |
| Sediment | Metamorph.40.46 | 0.1154241 | 0.001 | 0.078 | 0.008 |
| Sediment | Subadult | 0.1106084 | 0.001 | 0.078 | 0.008 |
| Sediment | Adult | 0.1625602 | 0.001 | 0.078 | 0.008 |
| Eggs.11.15 | Tadpole.27.29 | 0.0617868 | 0.001 | 0.078 | 0.008 |
| Eggs.11.15 | Metamorph.40.46 | 0.0969022 | 0.001 | 0.078 | 0.008 |
| Eggs.11.15 | Subadult | 0.1289319 | 0.001 | 0.078 | 0.008 |
| Eggs.11.15 | Adult | 0.1091131 | 0.001 | 0.078 | 0.008 |
| Tadpole.25.27 | Adult | 0.0902338 | 0.001 | 0.078 | 0.008 |
| Tadpole.27.29 | Adult | 0.0709375 | 0.001 | 0.078 | 0.008 |
| Metamorph.40.46 | Adult | 0.1078262 | 0.001 | 0.078 | 0.008 |
Q2: Are there transient versus symbiotic OTU’s on boreal toad skin?
Which fungi are more likely to be symbionts in that they are associated with the host as opposed to transiently moving through the host system?
To answer this question, I’m using ANCOM, an indicator species analysis for microbial sequence datasets. - Uses log-ratio transformation to somewhat account for compositionality of the data, which helps decrease the false discovery rate. - Uses K-W and wilcoxon tests to determine significance. - Seems to be more stringent than other methods, like LEfSe or macroecology indicator species analysis, in terms of how many OTU’s pass the thresholds. - Does not need rarefied data, but I am using my rarefied data here. - ANCOM is running a bunch of pairwise tests and the W stat is how many of those pairwise tests passed the “significance” threshold set in the initial script when I ran it.
Run ANCOM
# Ancom needs Sample.ID to be the header starting at A1, which means changing it in both files and transposing the taxa table.
# export_taxa_table(input = input_all,
# paste0(work_dir,"02_hypothesisbasedanalsysis/taxa_table_ancom.txt"),
# paste0(work_dir,"02_hypothesisbasedanalsysis/map_ancom.txt"))
# Did changes in excel, and changed OTU ID's to have the taxonomy and OTU # (for future reference)
# re-read in clean files
OTU_ancom <- read.delim(paste0(work_dir,"02_hypothesisbasedanalsysis/otu_ancom.txt"),
sep = "\t", header = T)
map_ancom <- read.delim(paste0(work_dir,"02_hypothesisbasedanalsysis/map_ancom.txt"),
sep = "\t", header = T)
# source the code, downloaded from: https://github.com/mech3132/MSC_code/blob/master/ANCOM_updated_code_MYC.R
source("/Users/alal3493/Documents/Reference_Info/ANCOM_MYC/ANCOM_updated_code_MYC_AA.R")- ANCOM run using sample type as toad or environmental samples
# run ancom with parameters for no repeated data, no adjustment (so it'll use a Wilcoxon test)
comparison_test_lifeenv <- ANCOM.main.myc(OTUdat = OTU_ancom,
Vardat = map_ancom,
adjusted = F,
repeated = F,
main.var = "Type",
adj.formula = NULL,
repeat.var = NULL,
longitudinal = F,
random.formula = NULL,
multcorr = 2,
sig = 0.05,
prev.cut = 0.90)
## [1] "Removing high zero samples..."
## [1] "Beginning logratio calculations..."
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] "Calculating critical statistic values..."
## [1] "Generating plots..."
## [1] "Done"
# This tells you if your data has a lot of zeroes/skew, as most microbial data do.
# everything below the red line gets removed.
comparison_test_lifeenv$PLot.zeroes
# This tells you which (colored) OTU's have a high effect size and significance
comparison_test_lifeenv$Plot.volcano
comparison_test_lifeenv$Plot.volcano.0.9
# this helps you determine if the critical w was chosen correctly. We want the valley between a bimodal distribution, in theory.
comparison_test_lifeenv$Plot.critical.w
# I think W=50 is good since it is very stringent, especially with the clumping of the points in the previous plots.
# ggsave(paste0(work_dir,"Figs_Tables/Wplot_lifeenv.png"),
# comparison_test_lifeenv$Plot.critical.w,
# dpi = 300, height = 5, width = 10)
# save(comparison_test_lifeenv,
# file = paste0(work_dir,"02_hypothesisbasedanalsysis/ANCOM_lifeenv.RData"),
# envir = .GlobalEnv)
# write.table(comparison_test_lifeenv$W.taxa, paste0(work_dir,"02_hypothesisbasedanalsysis/ANCOMlifeenv_taxabunds.txt"), sep="\t", quote=F)
# load(paste0(work_dir, "02_hypothesisbasedanalsysis/ANCOM_lifeenv.RData"))- ANCOM comparing each individual life stage to each environmental type.
list_env <- c("Sediment", "Water")
list_life <- c("Eggs","Tadpoles","Metamorphs","Subadults","Adults")
ALL_PAIRWISE <- list()
for (e in list_env) {
for (l in list_life) {
tryCatch({
# Mapping file filtered
mf_temp <- map_ancom %>%
filter(Life_Stage_Simplified %in% c(e, l))
# OTU file filtered
otu_temp <- input_all$data_loaded[ , which(colnames(input_all$data_loaded) %in% mf_temp[ , 1])]
# Transpose/ANCOM formatting
otu_temp_t <- as.data.frame(t(otu_temp)) %>%
rownames_to_column("Sample.ID")
mf_temp$Sample.ID
mf_temp$Type
colSums(otu_temp_t[,-1])
ALL_PAIRWISE[[paste0(e,"_",l)]] <- ANCOM.main.myc(OTUdat = otu_temp_t,
Vardat = mf_temp,
adjusted = F,
repeated = F,
main.var = "Type",
adj.formula = NULL,
repeat.var = NULL,
longitudinal = F,
random.formula = NULL,
multcorr = 2,
sig = 0.05,
prev.cut = 0.9)
}, error = function(e){cat("ERROR: ",conditionMessage(e), "\n")})
}
}
## [1] "Removing high zero samples..."
## [1] "Beginning logratio calculations..."
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] "Calculating critical statistic values..."
## [1] "Generating plots..."
## [1] "Done"
## [1] "Removing high zero samples..."
## [1] "Beginning logratio calculations..."
##
|
| | 0%ERROR: all observations are in the same group
## [1] "Removing high zero samples..."
## [1] "Beginning logratio calculations..."
##
|
| | 0%ERROR: all observations are in the same group
## [1] "Removing high zero samples..."
## [1] "Beginning logratio calculations..."
##
|
| | 0%ERROR: all observations are in the same group
## [1] "Removing high zero samples..."
## [1] "Beginning logratio calculations..."
##
|
| | 0%ERROR: all observations are in the same group
## [1] "Removing high zero samples..."
## [1] "Beginning logratio calculations..."
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%
## [1] "Calculating critical statistic values..."
## [1] "Generating plots..."
## [1] "Done"
## [1] "Removing high zero samples..."
## [1] "Beginning logratio calculations..."
##
|
| | 0%ERROR: all observations are in the same group
## [1] "Removing high zero samples..."
## [1] "Beginning logratio calculations..."
##
|
| | 0%ERROR: all observations are in the same group
## [1] "Removing high zero samples..."
## [1] "Beginning logratio calculations..."
##
|
| | 0%ERROR: all observations are in the same group
## [1] "Removing high zero samples..."
## [1] "Beginning logratio calculations..."
##
|
| | 0%ERROR: all observations are in the same group
# only sediment egg and water egg comparisons had significant differences
# the rest of the combos gave the message "all observations are in the same group"
# Calling sediment egg
ALL_PAIRWISE$Sediment_Eggs$PLot.zeroes
ALL_PAIRWISE$Sediment_Eggs$Plot.volcano.0.9
ALL_PAIRWISE$Sediment_Eggs$Plot.critical.w 
# really, could make it as low as 40 or 50, but let's stick to w = 124 ish as the cutoff to be extra stringent
# ggsave(paste0(work_dir,"Figs_Tables/Wplot_sed_egg.png"),
# ALL_PAIRWISE$Sediment_Eggs$Plot.critical.w,
# dpi = 300, height = 5, width = 10)
# calling water eggs
ALL_PAIRWISE$Water_Eggs$PLot.zeroes
ALL_PAIRWISE$Water_Eggs$Plot.volcano.0.9
ALL_PAIRWISE$Water_Eggs$Plot.critical.w 
# looks like W = 35 is best
# ggsave(paste0(work_dir,"Figs_Tables/Wplot_water_egg.png"),
# ALL_PAIRWISE$Water_Eggs$Plot.critical.w,
# dpi = 300, height = 5, width = 10)
# ran this three times and it gave the same results
# load(paste0(work_dir, "02_hypothesisbasedanalsysis/ANCOM_lifestage.RData"))
# save.image(file = paste0(work_dir, "02_hypothesisbasedanalsysis/FungalDevData_mainpub_workspace.RData"))Graph Volcano plot
- environment vs. toad samples
# load(file=paste0(work_dir,"02_hypothesisbasedanalsysis/ANCOM_lifeenv.RData"))
# make a volcano plot showing the associated microbes
W_frame <- comparison_test_lifeenv$W.taxa %>%
mutate(Enriched = case_when(W_stat < 50 | W_f < 10 ~ "Transient, not significant",
Environment_relAbund > Toad_relAbund ~ "Transient",
TRUE ~ "Symbiotic"))
# figure out order of ecol for legend
Ecol_list <- c("Symbiotic", "Transient", "Transient, not significant")
W_frame$Enriched <- factor(W_frame$Enriched, levels = Ecol_list)
# write.table(W_frame,
# paste0(work_dir,"02_hypothesisbasedanalsysis/indic_taxabunds_lifeenv.txt"),
# sep="\t", quote=F)
# graph volcano plot
volcano.plot <- ggplot(W_frame) +
geom_point(aes(x = W_f, y = W_stat, color = Enriched)) +
xlab("CLR F statistic (Effect size)") +
ylab("W statistic") +
labs(color = "Putative host relationship") +
scale_color_manual(values = c("#56B4E9", "#CC79A7", "#000000")) +
geom_hline(mapping = aes(yintercept = 50), col = "red", lty = 2) +
geom_vline(mapping = aes(xintercept = 10), col = "red", lty = 2)
volcano.plot
# ggsave(paste0(work_dir,"Figs_Tables/volcano_plot_color.png"),
# volcano.plot, dpi = 400, width = 7, height = 5)- Sediment vs. egg volcano plot
# make a volcano plot showing the associated microbes
W_frame_se <- ALL_PAIRWISE$Sediment_Eggs$W.taxa %>%
mutate(Enriched = case_when(W_stat < 123 | W_f < 10 ~ "Transient, not significant",
Environment_relAbund > Toad_relAbund ~ "Transient",
TRUE ~ "Symbiotic"))
# figure out order of ecol for legend
Ecol_list <- c("Symbiotic", "Transient", "Transient, not significant")
W_frame_se$Enriched <- factor(W_frame_se$Enriched, levels = Ecol_list)
# write.table(W_frame_se,
# paste0(work_dir,"02_hypothesisbasedanalsysis/indic_taxabunds_sedegg.txt"),
# sep = "\t", quote = F)
# graph volcano plot
volcano.plot.se <- ggplot(W_frame_se) +
geom_point(aes(x = W_f, y = W_stat, color = Enriched)) +
xlab("CLR F statistic (Effect size)") +
ylab("W statistic") +
labs(color = "Putative host relationship") +
scale_color_manual(values = c("#56B4E9", "#CC79A7", "#000000")) +
geom_hline(mapping = aes(yintercept = 123), col = "red", lty = 2) +
geom_vline(mapping = aes(xintercept = 10), col = "red", lty = 2)
volcano.plot.se
# ggsave(paste0(work_dir,"Figs_Tables/volcano_plot_sedegg.png"),
# volcano.plot.se, dpi = 400, width = 7, height = 5)- Water vs. egg volcano plot
# make a volcano plot showing the associated microbes
W_frame_we <- ALL_PAIRWISE$Water_Eggs$W.taxa %>%
mutate(Enriched = case_when(W_stat < 35 | W_f < 10 ~ "Transient, not significant",
Environment_relAbund > Toad_relAbund ~ "Transient",
TRUE ~ "Symbiotic"))
# figure out order of ecol for legend
Ecol_list <- c("Symbiotic", "Transient", "Transient, not significant")
W_frame_we$Enriched <- factor(W_frame_we$Enriched, levels = Ecol_list)
# write.table(W_frame_we,
# paste0(work_dir,"02_hypothesisbasedanalsysis/indic_taxabunds_wateregg.txt"),
# sep = "\t", quote = F)
# graph volcano plot
volcano.plot.we <- ggplot(W_frame_we) +
geom_point(aes(x = W_f, y = W_stat, color = Enriched)) +
xlab("CLR F statistic (Effect size)") +
ylab("W statistic") +
labs(color = "Putative host relationship") +
scale_color_manual(values = c("#56B4E9", "#000000")) +
geom_hline(mapping = aes(yintercept = 35), col = "red", lty = 2) +
geom_vline(mapping = aes(xintercept = 10), col = "red", lty = 2)
volcano.plot.we
# ggsave(paste0(work_dir,"Figs_Tables/volcano_plot_wateregg.png"),
# volcano.plot.we, dpi = 400, width = 7, height = 5)- create a three paneled version of all the above volcano plots
legend_volcano <- arrangeGrob(get_legend(volcano.plot))
volc_lifeenv_noleg <- volcano.plot + ggtitle("A") + theme(legend.position = "none")
volc_se_noleg <- volcano.plot.se + ggtitle("B") + theme(legend.position = "none")
volc_we_noleg <- volcano.plot.we + ggtitle("C") + theme(legend.position = "none")
# Put all into paneled graph with annotations and shared legend
volcano_fig <- ggpubr::ggarrange(volc_lifeenv_noleg,
volc_se_noleg,
volc_we_noleg,
legend_volcano,
ncol = 4, widths = c(5,5,5,2))
volcano_fig
# ggsave(paste0(work_dir,"Figs_Tables/volcano_combo.png"), volcano_fig,
# dpi = 500,
# width = 17, height = 3)A relative abundance dot plot using the data
# filter to only have ecologically relevant data
lifeenv_filt_ecol <- W_frame %>%
dplyr::filter(Enriched %in% c("Symbiotic", "Transient"))
se_filt_ecol <- W_frame_se %>%
dplyr::filter(Enriched %in% c("Symbiotic", "Transient"))
we_filt_ecol <- W_frame_we %>%
dplyr::filter(Enriched %in% c("Symbiotic", "Transient"))
merged_filt_ecol <- rbind(lifeenv_filt_ecol, se_filt_ecol, we_filt_ecol)
we_se_filt_ecol <- rbind(se_filt_ecol, we_filt_ecol)
# write.table(merged_filt_ecol,
# paste0(work_dir,"02_hypothesisbasedanalsysis/indic_taxabunds_combo.txt"),
# sep = "\t", quote = F)
rel_plot_lifeenv <- ggplot(lifeenv_filt_ecol) +
geom_point(aes(x = Environment_relAbund, y = Toad_relAbund, color = Enriched)) +
xlab("Relative Abundance on Environmental Samples") +
ylab("Relative Abundance on Toad Samples") +
labs(color = "Putative host relationship") +
xlim(0, 0.3) + ylim(0, 0.3) +
scale_color_manual(values = c("#56B4E9", "#CC79A7")) +
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed")
rel_plot_lifeenv
rel_plot_we_se <- ggplot(we_se_filt_ecol) +
geom_point(aes(x = Environment_relAbund, y = Toad_relAbund, color = Enriched)) +
xlab("Relative Abundance on Environmental Samples") +
ylab("Relative Abundance on Toad Samples") +
labs(color = "Putative host relationship") +
xlim(0, 0.075) + ylim(0, 0.075) +
scale_color_manual(values = c("#56B4E9", "#CC79A7")) +
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed")
rel_plot_we_se
# create paneled fig
legend_relplot <- arrangeGrob(get_legend(rel_plot_lifeenv))
relplot_le_noleg <- rel_plot_lifeenv + ggtitle("A") + theme(legend.position = "none")
relplot_wese_noleg <- rel_plot_we_se + ggtitle("B") + theme(legend.position = "none")
# Put all into paneled graph with annotations and shared legend
relplot_fig <- ggpubr::ggarrange(arrangeGrob(relplot_le_noleg,
relplot_wese_noleg),
legend_relplot,
ncol = 2, widths = c(4,3))
relplot_fig
# ggsave(paste0(work_dir,"Figs_Tables/relplot_fig.png"), relplot_fig,
# dpi = 300, width = 7, height = 7)Miscellaneous
Here, I’m making a graph showing where certain sample types fell out in the various processing steps.
drops <- read.delim(paste0(work_dir,"Figs_Tables/Metadata_dropouts_R.txt"))
#reorder the axes
all_sampletypes <- c("Upland_Soil", "Sediment", "Water", "Egg", "Tadpole", "Metamorph", "Subadult", "Adult", "Mouth", "Tadpole_Dead", "Glove", "Sterile_water", "NoTemplateControl")
drops$Life_Stage_Simplified <- factor(drops$Life_Stage_Simplified,
levels = all_sampletypes)
madetosteplist <- c("noDNA", "filtSeq", "filtR", "Final")
drops$Made_to_step <- factor(drops$Made_to_step, levels = madetosteplist)
#make graph
drop_graph <- ggplot(drops, aes(fill = Life_Stage_Simplified,
x = Made_to_step)) +
geom_bar(position = "dodge", stat = "count") +
labs(fill = "Sample Type", x = "Reasons for samples dropping out",
y = "Count") +
scale_x_discrete(labels = c("PCR did not amplify target gene",
"Filtered during sequence processing
(low reads or quality)",
"Filtered during R analysis
(some sample types not useful
for answering hypotheses)",
"Final dataset")) +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 14),
axis.title = element_text(size = 16),
axis.text.y = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_text(size = 16))
drop_graph
# ggsave("drop_graph.jpg", drop_graph, width=21, height=10)# save R workspace
# save.image(file = paste0(work_dir, "02_hypothesisbasedanalsysis/FungalDevData_mainpub_workspace.RData"))